Efficient Monte Carlo simulation of a glass forming binary mixture 
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We propose and use a novel, hybrid Monte Carlo algorithm that combines configurational bias 
particle swaps with parallel tempering. We use this new method to simulate a standard model of a 
glass forming binary mixture above and below the so-called mode-coupling temperature, Tmct- We 
find that an ansatz that was used previously to extrapolate thermodynamic quantities to tempera- 
tures below Tmct breaks down in the vicinity of the mode-coupling temperature. Thus, previous 
estimates of the so-called Kauzmann temperature need to be reexamined. Also, we find that the 
Adam-Gibbs relations D oc exp{—a/TSc) and r oc exp(&/rS'c), which connect the diffusion coeffi- 
cient D and the relaxation time r with the configurational entropy Sc, are valid for all temperatures 
for which the configurational and vibrational contributions to the free energy decouple. 
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Understanding the nature of the glass transition has 
been of great interest for several decades. One of the 
earliestparadigms P, Q , that has recently been reformu- 
lated y and subsequently received considerable atten- 
tion, assumes the existence of an "ideal" glass transition 
which occurs when the entropy of the supercooled liq- 
uid becomes equal to the entropy of a disordered solid. 
Thisparadigm has stimulated several simulational stud- 
ies |5| which have confirmed both the qualitative de- 
scription of supercooled liquids' dynamics and even the 
quantitative, numerical prediction of the transition tem- 
perature, the so-called Kauzmann temperature Tk, for a 
popular model of a glass forming binary mixture. How- 
ever, these simulational investigations were restricted to 
temperatures above the so-called mode-coupling temper- 
ature, Tmct- Thus, in these previous studies, in or- 
der to investigate the equality of the liquid and disor- 
dered solid entropies, one had to extrapolate higher tem- 
perature data obtained directly from simulations to sig- 
nificantly lower temperatures (for example, for the bi- 
nary mixture studied in Refs. |3j it was found that 
Tmct /Tk ~ 1.45). The commonly used extrapolations 
were questioned in a recent investigation of the configu- 
rational entropy |^. This study is unique in that it ac- 
cessed directly the low temperature region below Tmct- 
It showed that the commonly used extrapolations break 
down below Tmct and the existence of the Kauzmann 
temperature was put in doubt. Reference used a novel 
density-of-states Monte Carlo method. Also, it used a bi- 
nary mixture model that was somewhat different 0| than 
the model used in most of the previous simulations. Fi- 
nally, although several system sizes were considered in 
Ref. 1^, the largest system (216 particles) was signifi- 
cantly smaller than systems used in previous studies (the 
largest ones had 1000 particles). 

The goal of our investigation was to address the ques- 
tion of the existence of the Kauzmann temperature for 
the original model studied in Refs. 0, 0|. Since the 
long relaxation times below Tmct makes equilibration 



of molecular and Brownian dynamics simulations very 
difhcult, we propose and use a novel, specialized Monte 
Carlo algorithm designed to decrease the time between 
independent measurements of strongly supercooled liq- 
uids. This new method allows us to obtain accurate ther- 
modynamic quantities for temperatures below the mode- 
coupling temperature. We find that the previously used 
extrapolations of thermodynamic quantities to estimate 
the Kauzmann temperature break down in the vicinity of 
the mode-coupling temperature. Furthermore, by using 
the results of Brownian dynamics simulations, we demon- 
strate that the Adam-Gibbs relations are valid for low 
temperatures. We propose using the Adam-Gibbs rela- 
tion and the results of the Monte Carlo simulation to 
predict the diffusion coefRcient. 

To directly access temperatures below the mode- 
coupling temperature, we combined non-local trial moves 
with parallel tempering. Non-local trial moves, identity 
exchanges, have been effective in speeding up simulations 
with different size particles However, high den- 

sity and large size disparity drastically reduce the accep- 
tance rate of identity exchanges , which decreases their 
usefulness in dense glass forming liquids at low temper- 
atures. To overcome this, we used parallel tempering, 
which attempts to exchange replicas of the system which 
are being simulated at different temperatures. This al- 
lows configurations which are trapped in a deep potential 
energy minimum at low temperatures to change consider- 
ably by being simulated at higher temperatures. A vari- 
ation of this technique, which combines replica exchange 
with molecular dynamics, has been shown to spee d up 
the equilibration of the system studied in this work [lO| . 

We simulated an 80:20 binary mixture introduced by 
W. Kob and H. Andersen The interaction poten- 

tial is Vafsinj) = 4ea/3 [{(Tafs/rij)^'^ - (aap/rij)^] where 
a,/3 G {A,B}, eAA = 1-0, e^s = 1-5, ess = 0.5, 
aAA — 1-0, (JAB — 0.8, and ctbb — 0.88. The results 
are presented in reduced units with eAA and uaa being 
the units of energy and length, respectively. We simu- 
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lated 1000 particles in a fixed cubic box with a box length 
of 9.4. We performed a parallel tempering Monte Carlo 
simulation and a series of Brownian dynamics simula- 
tions. The Monte Carlo simulation was performed using 
the following set of temperatures, T — 0.62, 0.59, 0.56, 
0.53, 0.50, 0.48, 0.46, 0.44, 0.43, 0.42, 0.41, and 0.40. 
The Brownian dynamics simulations were performed at 
T = 5.0, 3.0, 2.0, 1.5, 1.0, 0.9, 0.8, 0.6, 0.55, 0.5, 0.47, 
0.45, and 0.44. 

The details and results of the Brownian dyna mics sim- 
ulation have been presented elsewhere 0, 0| . Here we 
briefly describe the Monte Carlo simulation. We uti- 
lize three trial moves: a standard, local single particle 
displacement, a configurational bias particle swap, and 
parallel tempering. The configuration bias particle ex- 
change attempts to swap two particles of different sizes. 
The smaller particle will always fit in the space left by 
the larger particle, but, since the density is high, the 
converse is rarely true. To increase the acceptance rate, 
50 trial configurations are explored for the larger particle 
around the former position of the smaller particle. One of 
the trial positions is chosen with a probability which de- 
pends on the potential energy, and the move is accepted 
with a probability so that detailed balance is maintained 
The configurational bias increases the acceptance 
rate, but the acceptance rate for the particle swaps is still 
very small, around 10~^ at T = 0.5, which is the lowest 
temperature in which the swaps were attempted. It has 
been shown that identity exchange can decrease the equi- 
libration time of a simulation dramatically even if 
the acceptance rate is very small. Parallel tempering 
consists of an attempted exchange of particle positions 
between adjacent temperatures Thus, the configu- 
rational bias particle swaps can decrease the correlation 
time for the configurations which begin the simulation at 
the temperatures in which the swaps are not attempted. 

To estimate the efficiency of our algorithm, we com- 
pared the energy correlation time measured in Monte 
Carlo moves per particle with the energy correlation time 
in the Brownian Dynamics simulation measured in Brow- 
nian Dynamics time steps. We found that the former 
time increases slower with decreasing temperature than 
the latter time. The energy correlation time for the 
Monte Carlo simulation at T = 0.4 corresponds to gener- 
ating one statistically independent configuration approx- 
imately every five days for a parallel algorithm using four 
threads on a hyperthreaded dual- processor 3.2 GHz Pen- 
tium workstation. 

We used a variety of different checks for equilibration: 
monitoring the running average of the potential energy 
U{N) = {1/N) XliLi Ui where Ui is the potential energy 
for step i {U{N) does not show any systematic drift), 
comparing specific heat calculated using the derivative 
of the energy and energy fluctuations (they agree, see 
Fig. |2Jl, comparing the temperature assumed in the sim- 
ulation algorithm with the so-called configurational tem- 
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FIG. 1: The re- weighted potential energy proba- 
bility distributions. The solid line is a Gaussian 
distribution ^l/{2TvkBCv,pot{T)T^) GXp ( — ( t/ pot — 

(Upot))^ /{2kBCv,pot(T)T^) for T = 0.44 and T = 0.50. 



perature (they agree), etc. Here we discuss 

in some detail a stringent equilibration test introduced 
by Yamamoto and Kob 10]. This so-called re- weighting 
procedure relies upon the fact that an equilibrium dis- 
tribution of the energy for a temperature Ti can be de- 
termined from an equilibrated simulation at a different 
temperature Tq. Shown in Fig. ^ and Fig. ^ are the 
re- weighted probability distributions of the potential en- 
ergy for T — 0.44, and 0.5, respectively. The filled circles 
in the figures are the probability distributions calculated 
from simulation runs at T = 0.44 and 0.5, respectively. 
The other symbols are the re-weighted probability dis- 
tributions, and the solid lines are Gaussian distributions 
determined from the average energies and the specific 
heats calculated at the respective temperatures. Note 
that there is a slight systematic deviation from a Gaus- 
sian distribution at the lowest temperatures. This devi- 
ation is within the uncertainty of the calculation, thus it 
is not clear whether it has any significance. In any case, 
it should be noted that the re-weighted probability dis- 
tributions still superimpose very well. Thus, at all tem- 
peratures there is good overlap between the re-weighted 
distributions and this fact provides strong evidence that 
the Monte Carlo simulation has properly equilibrated. 

The specific heat calculated from energy fluctuations 
and from the derivative of the energy determined from 
the Brownian dynamics and the Monte Carlo simulation 
is shown in Fig. El Note that a correction due to the 
finite simulation time [tR 0| has been applied to 
the specific heat calculated from the energy fluctuations. 
The agreement between the Monte Carlo, the Brownian 
dynamics, and the two methods of calculating the spe- 
cific heat is very good except for the derivative of the 
energy for the Brownian dynamics simulation at the low- 
est temperature. There is a peak in the specific heat 
around T « 0.45, which is close to the usually cited 
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FIG. 2: The specific lieat per particle as a function of T~^^^ 
calculated from energy fluctuations (closed symbols) and the 
derivative of the energy (open symbols) for the Brownian dy- 
namics (squares) and the Monte Carlo (circles). The solid line 
is Cv/N = 0.6aT~^^^ + 1.5 where a was obtained from fitting 
the average potential energy (Upot) to a function of the form 
aT^^^ _l_ I) (^^^jjg fl, ajjfj 5 parameters obtained from the fit are 
a — —8.6547 and b — 2.6362). Inset: the average potential 
energy per particle as a function of T'^''^; the solid line is a 
aT^/^ + b fit. 



mode-coupling temperature Tmct ~ 0.435 |ll| deter- 
mined from fits to the diffusion coefficient and the relax- 
ation time. It is clear from FigJ^hat the extrapolation 
that was used in prior studies 0, |a] is violated below the 
mode-coupling temperature. A similar violation was ob- 
served in a recent investigation of a different binary mix- 
ture 0. Recall that the former studies only simulated 
systems at temperatures higher than T^cTi whereas the 
latter one was able to access temperatures below Tmct- 
In Fig. 13^ we show the total liquid entropy, S and 
the so-called disordered solid entropy, Smb- To find the 
total entropy S we performed a standard thermodynamic 
integration along the T = 5.0 isotherm and then along 
the Vq = (9.4)3 isochore. For T > 0.62 we utilized a 
commonly used fit for the specific heat shown as a solid 
line in Fig. 13 For T < 0.62 wc numerically integrated 
the specific heat. The results of the numerical integration 
are shown as triangles in Fig. |3Ji. The solid line going 
through the triangles for T > 0.45 and slightly deviating 
from them for T < 0.45 is the standard extrapolation 
0, m that relies on using the fit shown as a solid line in 
Fig. 121 for all temperatures. 

To evaluate the disordered solid entropy we followed 
the procedure used in previous studies 0, 0, |U 
E3 . E^ . First, we determined the inherent structures 
by quenching 500-1000 configurations. Then we checked 
that for r < 0.62 the system can be described by a con- 
figurational and a vibrational part Q . We calculated the 
entropy of the disordered harmonic solid by diagonalizing 
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FIG. 3: (a) The total entropy per particle S/N (triangles) 
and the disordered solid entropy per particle Svit/N. The 
latter entropy has been calculated using the inherent struc- 
tures obtained from the Monte Carlo (circles) and the Brow- 
nian dynamics (squares) simulations. The solid lines are the 
commonly used extrapolations described in the text, (b) The 
configurational entropy per particle, Sc/N — S/N — Svib/N. 
The solid line is the extrapolation obtained from the diflerence 
of the extrapolations shown in (a) . The error bars for S/N, 
Svib/N, and Sc/N are smaller than the size of the symbols. 



the Hessian matrix calculated at the inherent structures 
and determined the vibrational frequencies uji. Next, 
we calculated the vibrational contribution to the entropy 

Svib = '^[1 ^ ln(/3ftwi)]^ where < • >' denotes an 

average over the inherent structures (note that through- 
out this paper we set Planck's constant equal to one). 
The results are shown as circles and squares in Fig. 
The solid line going through the circles and squares is 

obtained by fitting ^X]^=5i~'^ l^(^j)]^ ^o a polynomial in 
T of degree 2. This quantity, which is the contribution to 
Syib that originates from the the vibrational frequencies, 
is almost temperature-independent. 

In Fig. we show the configurational entropy Sc = 
S — Smb- Note that Sc was only calculated for temper- 
atures for which the the system can be divided into a 
configurational and a vibrational part, i.e. for T < 0.62. 
In Fig. IJh we also compare our results with the previ- 
ously used extrapolation of the configurational entropy. 
This extrapolation results in the Kauzmann temperature 
TV = 0.29, which is very close to previous estimates 
[J, y ■ Note that below the mode-couphng temperature 
our results deviate from this extrapolation. Thus, the 
previous estimates of the Kauzmann temperature have 
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FIG. 4: Test of the Adam Gibbs relations D oc exp{—a/TSc) 
(solid symbols and right vertical axis) and r oc exp{b/TSc) 
(open symbols and left vertical axis) for the A (circles) and 
the B particles (squares). Inset: The diffusion coefficient de- 
termined directly from the Brownian dynamics simulations 
(closed symbols) and predicted using the Adam-Gibbs re- 
lation (open symbols) for the A (circles) and B particles 
(squares). The Adam-Gibbs prediction uses Sc determined 
from the Monte Carlo simulations. The dashed curves are 
mode-coupling theory power law fits a{T — 0.435)'' over the 
temperature range 0.8 < T < 0.5. 

to be reexamined. The peak in the specific heat indi- 
cates that the total entropy is larger than the previous 
estimate, and that the Kauzmann temperature, if it ex- 
ists at all, is lower. Our results are consistent with those 
obtained in Ref. for a different binary mixture. 

We examined the Adam-Gibbs relations D cx 
exp{-a/TSc) and r oc exp(6/rS'c) l3- Shown in Fig. H 
are the logarithms of the diffusion coefficients and the 
relaxation times (the relaxation times are obtained from 
the self-intermediate scattering function in the usual way 
mill) plotted as a function of 1/TSc for T < 0.6. The 
diffusion coefficients and the relaxation times are deter- 
mined from the Brownian dynamics simulations. The 
configurational entropy is determined from the Monte 
Carlo simulation and the Brownian dynamics simula- 
tions. Although some curvature can be seen in the re- 
laxation time data, the fits are very good and verify the 
Adams-Gibbs relation for this temperature range. As- 
suming that the Adam-Gibbs relation holds at lower tem- 
peratures, we can predict the diffusion coefficient down to 
T = 0.4 from the results of the Monte Carlo simulation; 
see the inset in Fig. ^ 

In summary, we propose and use a novel, hybrid 
Monte Carlo algorithm which allows accurate calcula- 
tion of equilibrium quantities below the mode-coupling 
temperature. Our method combines non-local configu- 
ration bias particle swaps and parallel tempering. The 
results obtained with this algorithm put in doubt the 
commonly used extrapolation of the supercooled liquid 



entropy and previous estimates of the Kauzmann tem- 
perature. Moreover, we show that the Adam-Gibbs re- 
lations D (X exp{~a/TSc) and t oc exp{b/TSc) hold for 
all temperatures for which the configurational and vibra- 
tional contributions to the free energy decouple. 
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